Den här inlämningsuppgiften förutsätter att följande paket finns installerade:
mosaic
dplyr
geosphere
leaflet
Paket kan installeras via kommandot install.packages('packagename'), där 'packagename' är namnet på paketet, t.ex 'mosaic'.
Introduktion
I den första inlämningsuppgiften ska ni självständigt i grupper om tre analysera ett dataset i programmeringsspråket R. Till skillnad från datorlaborationerna finns det minimalt med kodexempel. Datorlaborationerna går igenom de flesta momenten som behandlas i inlämningsuppgiften, så se till att göra klart dessa innan.
Instruktioner
I denna inlämningsuppgift ska ni analysera ett datamaterial som beskriver ca 500 distrikt i Boston år 1970. Datasetet är en modifierad version1 av originaldata2 som användes i en studie3 där författarna predikterade medianhuspriser i olika distrikt med hjälp av en rad förklaringsvariabler.
Följande variabler finns i datasetet boston_census_data.Rdata (ladda ner), som innehåller 480 observationer. Notera att en observation motsvarar ett distrikt:
zoned_25k_p: Andel av stadsdelens bostadsmark ämnad för marklotter större än 25000 kvadratfot.
indust_p: Andel tunnland ägd av företag utanför detaljhandel.
borders_charles: Charles River dummy variabel (= 1 om området gränsar till floden, 0 annars).
NOx: Koncentration av kväveoxider (andelar per 10 miljon).
n_rooms_avg: Genomsnitt antal rum i ägda bostäder.
before_1940_p: Andel ägda bostäder byggda före 1940.
employ_dist: Viktat avstånd till fem arbetsförmedlingscentra i Boston.
radial_access: Index som mäter tillgång till stadsmotorvägar.
tax_rate: Fastighetsskatt per 10000 USD.
pupil_teacher_ratio: Lärartäthet mätt som elev per lärare.
lower_stat_pct: Procentandel med låg socioekonomisk status i termer av utbildning eller arbete.
dist_fenway_park: Avstånd till stadion Fenway Park.
Inlämningsuppgiften ska lämnas in i form av ett html-dokument genererat av Quarto. Kontrollera att ni inte får några felmeddelande när du genererar HTML-dokumentet. Läs sedan igenom HTML-dokumentet noggrant innan ni lämnar in det. Använd tydliga figurer och namnge axlarna med tydliga variabelnamn. Glöm inte att skriva era namn i början av dokumentet, där det nu står Namn 1, Namn 2 och Namn 3.
0. Ladda in data
💪 Uppgift 0.1
Ladda in datasetet Boston_census_data med följande kod.
Uppgift 0.1 - Svar
# Write your code hereload(file =url("https://github.com/StatisticsSU/SDA1/blob/main/assignments/assignment1/Boston_census_data.RData?raw=true"))
1. Kriminalitet i Boston
I detta avsnitt ska ni analysera kriminaliteten i Boston med hjälp av variabeln crime_rate.
💪 Uppgift 1.1
Vad kan man generellt säga om kriminaliteten i Boston 1970? Använd lämpliga figurer och mått för att ge en beskrivning.
Uppgift 1.1 - Svar
Variabeln crime_rate mäter brott per 1000 invånare. Vi börjar med att använda kommandot “favstats” för att få en bra överblick över variabeln. Favstats säger oss att det minsta värdet är 0.00632 och 88.9762 är det högsta, det ger en varians på 88,96988 vilket i sammanhanget är väldigt mycket och säger oss att brottsligheten skiljer sig väldigt mycket åt mellan olika distrikt. Medelvärdet är 3.66526 med en standardavvikelse på 8.746 medan Medianen enbart är 0.253715 med en IQR på 3.599. I och med att datasetet är snedvridet så anser vi att medianen är ett bättre mått för denna variabel.
Genom att använda en boxplot bekräftar det tanken på att datasetet är väldigt snedvridet då vi ser väldigt många outliers. Vi kan även se om vi gör ett histogram istället att majoriteten av distrikten, över 80%, har en crime_rate på mellan 0 och 10 brott per tusen invånare, vilket med andra ord betyder att de flesta distrikten har en relativt låg crime_rate medan några få distrikt har en väldigt hög crime_rate.
# Write your code heresuppressMessages(library(mosaic))suppressMessages(library(dplyr))suppressMessages(library(geosphere))suppressMessages(library(leaflet))favstats(~crime_rate, data = Boston_census_data)
min Q1 median Q3 max mean sd n missing
0.00632 0.0829725 0.253715 3.681942 88.9762 3.66526 8.746156 480 0
bwplot(~crime_rate, data=Boston_census_data)
histogram(~crime_rate, data = Boston_census_data, breaks =10)
💪 Uppgift 1.2
Distrikten tillhör olika stadsdelar som anges i den kategoriska variabeln town? Det finns 88 olika sådana stadsdelar (towns).
Skiljer sig brottsligheten åt mellan de olika stadsdelarna? Undersök stadsdelarna Boston East Boston, Boston Downtown, Cambridge, samt ytterligare en stadsdel som ni själva väljer. Besvara frågan med hjälp av lämpligt valda figurer och statistiska mått.
Tip
Innan ni påbörjar analysen, skapa en ny data frame som enbart innehåller de stadsdelar som ni vill jämföra. Det kan göras till exempel med funktionen filter().
Uppgift 1.2 - Svar
Efter att vi har skapat vårt filter som enbart innehåller de utvalda stadsdelarna, n=63 istället för 480, kan precis som på förra deluppgiften börja med kommandot favstats. Där ser vi att variansen fortfarande är väldigt stor 51.1358 - 0.29819 = 50.83761 samt att medianen och medelvärdet skiljer sig åt, median = 2.3139 och medelvärdet = 6.261704. Det säger oss att fördelningen fortfarande är skev, både medelvärdet och medianen är högre här än om vi tittar på alla observationer.
Tittar vi på histogrammet ser vi tydligt att den är högerskev med en markant outlier med en crime_rate på 50 per 1000 invånare. Fördelningen är unimodal med ett typvärde på en crime_rate mellan 0 och 10 per 1000 invånare.
I boxploten ser vi att crime_raten är minst i Newton med väldigt liten spridning medan den störst i Boston Downtown med en median på runt 17-18, det finns också en outliner på över 50.
# Write your code hereBoston_census_new <-filter(Boston_census_data, town =="Boston East Boston"|town =="Boston Downtown"|town =="Cambridge"|town =="Newton")favstats(~crime_rate, data = Boston_census_new)
min Q1 median Q3 max mean sd n missing
0.29819 0.61913 2.3139 8.07211 51.1358 6.261704 9.129887 63 0
histogram(~crime_rate, data = Boston_census_new)
bwplot(crime_rate ~ town, data = Boston_census_new)
💪 Uppgift 1.3
Vilka två variabler i datasetet Boston_census_data korrelerar mest med brottslighet? Beskriv sambandet mellan brottslighet och var och en av dessa två variabler.
Tip
Kom ihåg att korrelation mäter det linjära sambandet mellan numeriska variabler.
Uppgift 1.3 - Svar
Vi skapar ett nytt dataset som enbart innehåller numeriska variabler som vi döper till “Boston_census_data_15_variables”. Sen för att hitta vilka två variabler som korrelerar mest med brottslighet gör vi en correlation matrix, denna är dock svåröversiktlig så då är det bättre att göra en corrplot, där man tydligare ser vad som korrelerar med vad. Detta kan man senare dubbelkolla med correlation matrixen, ifall vissa pluppar är snarlika varandra.
Vi hittar 3 variabler som ser ut att korrelera mest med crime_rate: radial_access, median_home_value och tax_rate. Vi använder funktionen cor() för att se vilka två som har högst samband. radial_access: 0.62, median_home_value: -0.45 och tax_rate: 0.58.
Radial_access och crime_rate har ett relativt starkt positivt samband på 0.62, det innebär ju större tillgång till stadsmotorvägar desto högre crime_rate.
tax_rate och crime_rate har också ett postivt samband på 0.58, det innebär att ju högre tax_rate desto mer brott begås per 1000 invånare.
# Write your code hereBoston_census_data_15_variables <- Boston_census_data[, c("longitude","latitude","median_home_value","crime_rate","zoned_25k_p","indust_p","NOx","n_rooms_avg","before_1940_p","employ_dist","radial_access","tax_rate","pupil_teacher_ratio","lower_stat_pct","dist_fenway_park")]correlation_matrix_Boston <-cor(Boston_census_data_15_variables)round(correlation_matrix_Boston, 3)
cor(radial_access ~ crime_rate, data = Boston_census_data)
[1] 0.623462
cor(median_home_value ~ crime_rate, data = Boston_census_data)
[1] -0.4496238
cor(tax_rate ~ crime_rate, data = Boston_census_data)
[1] 0.5802444
2. Fastighetsskatt i Boston
I detta avsnitt ska ni analysera fastighetsskatten i Boston med hjälp av variabeln tax_rate.
💪 Uppgift 2.1
Vad kan man generellt säga om fastighetsskatten i distrikten? Använd lämpliga figurer och mått för att beskriva fördelningen.
Uppgift 2.1 - Svar
Fastighetsskatten varierar mellan 187 och 711 per 1000 USD, vilket ger en varians på 524. Medianen och Medelvärdet skiljer sig inte lika mycket som det gjorde på brottsligheten, utan ligger här på 330 respektive 409.33. Standardavvikelsen är 168.5777 och IQR är 385.
Om vi gör ett histogram ser vi att det är en bimodal fördelning samt att det ser ut att vara två olika grupper i och med att det är ett glapp mellan staplarna. Typvärderna ligger runt 300 och 650. Antagligen finns det ett samband med en annan variabel som avgör om man betalar lägre eller högre fastighetsskatt.
# Write your code herefavstats(~tax_rate, data=Boston_census_data)
min Q1 median Q3 max mean sd n missing
187 280.75 330 666 711 409.3271 168.5777 480 0
histogram(~ tax_rate, data = Boston_census_data, col ="navyblue", type="count")
💪 Uppgift 2.2
Låt oss skapa en ny variabel cat_tax som anger om ett distrikt betalar låg (low), medel (medium), eller hög (high) fastighetsskatt. Vi definerar skattekategorierna enligt
low: tax_rate\(\leq\) 250,
medium: 250 \(<\)tax_rate\(\leq\) 400,
high: tax_rate\(>\) 400.
Följande kod skapar och lägger till variabeln cat_tax i Boston_census_data
Finns det ett samband mellan vilken skattekategori ett distrikt tillhör och dess angränsning till Charles River? Förklara med hjälp av lämplig tabell och figur.
Uppgift 2.2 - Svar
Som man ser i stapeldiagrammet skiljer sig skattenivåerna inte så mycket åt mellan om man gränsar till Charles River eller inte. Det man kan se är att det är mer troligt att man har men medium-tax rate om man gränsar till Charles River än om man inte gör det.
Fördelningen mellan skattenivåer om man inte gränsar till Charles River eller om man gör det kan vi se i korstabellen. Det är större sannolikhet att man inte gränsar till Charles River om man har låg (13,3% mot 10,3%) eller hög (39,69% mot 34,48%) skattekategori.
# Write your code heretally(cat_tax ~ borders_charles, data = Boston_census_data, margins=TRUE, format="percent")
borders_charles
cat_tax 0 1
Low 13.30377 10.34483
Medium 47.00665 55.17241
High 39.68958 34.48276
Total 100.00000 100.00000
bargraph(~cat_tax, groups=borders_charles, data=Boston_census_data, type="percent", main ="Skattekategorier baserat på angräsning till Charles River")
💪 Uppgift 2.3
Hur många procent av alla distrikt i vår data ligger i angränsning till Charles River och tillhör en hög skattekategori? Hur stor andel av distrikten med hög skatt ligger inte i angränsning till Charles River?
Uppgift 2.3 - Svar
Nedan ser vi först en marginalfördelningstabell. Där kan vi se att 2,08% av alla distrikt i vår data angränsar till River Charles (1) och tillhör en hög skattekategori (High).
Vi ser sedan en betingad fördelningstabell, som är uppdelad på skattekategorierna. I den kan vi se att av de distrikt som har en hög skattekategori är det 94,71% som inte ligger i angränsning till River Charles (0).
borders_charles
cat_tax 0 1 Total
Low 12.500000 0.625000 13.125000
Medium 44.166667 3.333333 47.500000
High 37.291667 2.083333 39.375000
Total 93.958333 6.041667 100.000000
tally(borders_charles ~ cat_tax, data = Boston_census_data, margins=TRUE, format="percent")
cat_tax
borders_charles Low Medium High
0 95.238095 92.982456 94.708995
1 4.761905 7.017544 5.291005
Total 100.000000 100.000000 100.000000
💪 Uppgift 2.4
Vilka två variabler i datasetet Boston_census_data korrelerar starkast med tax_rate? Beskriv det parvisa sambandet mellan tax_rate och var och en av dessa två variabler. Vad kan vi säga om kausalitet för vart och ett av sambanden?
Tip
Kom ihåg att korrelation är ett mått på linjära samband mellan numeriska variabler.
Uppgift 2.4 - Svar
Skriv svaret här.
Vi skapar ett nytt dataset som enbart innehåller numeriska variabler som vi döper till “Boston_census_data_15_variables”. Sen för att hitta vilka två variabler som korrelerar mest med fastighetsskatt gör vi en correlation matrix, denna är dock svåröversiktlig så då är det bättre att göra en corrplot, där man tydligare ser vad som korrelerar med vad. Detta kan man senare dubbelkolla med correlation matrixen, ifall vissa pluppar är snarlika varandra.
Man ser ganska tydligt att det är två variabler som verkar korrelera mest med tax_rate och det är radial_access och indust_p. Vi använder cor()-funktionen för att se korrelationerna.
Radial_access och tax_rate har ett väldigt starkt positivt samband på 0.91. Det verkar finnas ett tydligt samband mellan större tillgång till stadsmotorvägar och distriktens fastighetsskatt.
Även indust_p och tax_rate har ett starkt positivt samband på 0.72. indust_p står för andel tunnland ägd av företag utanför detaljhandeln. Det kan man tolka att ju mer tunnland som ägs av företag utanför detaljhandeln, desto högre tax_rate.
Vi kan inte säga att någon av dessa är kausala samband, vi vet inte med säkerhet att det är t ex tillgången till stadsmotorvägar som påverkar fastighetsskatten, utan det kan alltid finnas en dold variabel utanför vårt dataset som påverkar både, och bara får det att se ut som att det finns ett samband.
# Write your code hereBoston_census_data_15_variables <- Boston_census_data[, c("longitude","latitude","median_home_value","crime_rate","zoned_25k_p","indust_p","NOx","n_rooms_avg","before_1940_p","employ_dist","radial_access","tax_rate","pupil_teacher_ratio","lower_stat_pct","dist_fenway_park")]suppressMessages(library(corrplot))corrplot(correlation_matrix_Boston)
cor(radial_access ~ tax_rate, data = Boston_census_data)
[1] 0.9099416
cor(indust_p ~ tax_rate, data = Boston_census_data)
[1] 0.7171041
3. Avstånd till Fenway park
I detta avsnitt ska ni undersöka variabeln dist_fenway_park, som mäter avståndet mellan ett distrikt och Fenway park (stadion där basebollslaget Boston Red Sox spelar sina hemmamatcher).
Vi kan visualisera Fenway park och distrikten på en karta med hjälp av R-paketet leaflet. Följande kod visar platsen för Fenway park och distrikten för observationerna 30 och 45.
library(leaflet) # Install if not availablefenway_park_lat_long <-c(42.346462, -71.097250) # latitude and longitude for Fenway_parkBoston_map <-leaflet() %>%addTiles() %>%addMarkers(lat = fenway_park_lat_long[1], lng = fenway_park_lat_long[2], popup="Fenway park") %>%addMarkers(lat = Boston_census_data$latitude[30], lng = Boston_census_data$longitude[30], popup="Observation 30") %>%addMarkers(lat = Boston_census_data$latitude[45], lng = Boston_census_data$longitude[45], popup="Observation 45") Boston_map # Show interactive map
💪 Uppgift 3.1
Vilket distrikt i vår data har längst respektive kortast avstånd till Fenway park? Markera ut dessa distrikt i en interaktiv karta tillsammans med Fenway park.
Uppgift 3.1 - Svar
Det första vi måste göra är att få fram vilka distrikt som har längst respektive kortast avstånd till Fenway park. Det får vi fram genom att sortera datan på “dist_fenway_park”-variablen. Om vi använder funktionen head() kan vi se distrikten på först raden. Wilmington har kortast avstånd på 887.9 medan Marshfield har längst på 33638.4.
För att få ut dessa på kartan, betyder vi ut till våra nya variabler för den sorterade datan och väljer den första raden i respektive.
Finns det ett samband mellan dist_fenway_park och crime_rate?
Uppgift 3.2 - Svar
Det finns ett svagt negativt samband mellan avståndet till Fenway park och crime_rate. Korrelationen antyder att närmare Fenway park (lägre avstånd) desto mer brottslighet (högre crime_rate). Dock är sambandet svagt så det kan likaväl bero på slumpen.
# Write your code herecor(dist_fenway_park ~ crime_rate, data = Boston_census_data)
[1] -0.1228651
4. Enkel linjär regression
I detta avsnitt ska ni anpassa och tolka enkla linjära regressionsmodeller.
💪 Uppgift 4.1
Anpassa en linjär regression med responsvariabeln NOx och den förklarande variabeln employ_dist. Rita den anpassade regressionslinjen tillsammans med data i en lämplig figur. Beskriv resultaten och tolka modellen. Utför en modellvalidering via en residualanalys och kommentera modellens lämplighet. Om modellen inte anses lämplig, vilka antaganden har inte varit uppfyllda?
Uppgift 4.1 - Svar
Vår responsvariabel,y, är NOx, vilket är koncentration av kväveoxider (andelar per 10 miljoner). Variabeln som vi ska använda för att försöka förklara variantionen av NOx är employ_dist, x, vilket är det viktade avståndet till 5 arbetsförmedlingscentra i Boston.
Utifrån regressionsmodellen kan vi se att interceptet är 0.7176, vilket är det värde som NOx har när employ_dist är 0. Lutningen, slopen är -0.0426 vilket behöver att varje gång employ_dist ökar med en enhet så minskar NOx med -0.0426. Detta kan man tolka som att ju längre avståndet är till arbetsförmedlingscentra desto mindre är koncentrationen av kväveoxider.
R-squared är 0.5861, vilket betyder att employ_dist kan förklara ca 59% av variationen i NOx.
Modellvalidering
Utifrån modellvalidering kan vi se att modellen inte uppfyller något av kriterierna. I residualplotten kan vi se att residulernas varians inte är konstant, dem är koncentrerade på den högra sidan av modellen samt har några outliners uppe i det vänstra hörnet.
Den är inte heller helt normalfördelad, då den inte följer det röda strecket i början eller slutet, även om de gör det på det stora hela.
# Write your code hereregressionmodel_NOx_employdist <-lm(NOx ~ employ_dist, data = Boston_census_data)summary_model <-summary(regressionmodel_NOx_employdist)summary_model
Call:
lm(formula = NOx ~ employ_dist, data = Boston_census_data)
Residuals:
Min 1Q Median 3Q Max
-0.12358 -0.05352 -0.01220 0.04348 0.22868
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.717614 0.007109 100.95 <2e-16 ***
employ_dist -0.042641 0.001639 -26.02 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07526 on 478 degrees of freedom
Multiple R-squared: 0.5861, Adjusted R-squared: 0.5852
F-statistic: 676.9 on 1 and 478 DF, p-value: < 2.2e-16
plot(NOx ~ employ_dist, data = Boston_census_data, col ="cornflowerblue", ylim =c(0, 1))y_hat <-predict(regressionmodel_NOx_employdist)lines(Boston_census_data$employ_dist, y_hat, type ="p", col ="lightcoral")abline(regressionmodel_NOx_employdist, col ="lightcoral")legend(x ="topleft", pch =c(1, 1, NA), lty =c(NA, NA, 1), col =c("cornflowerblue", "lightcoral", "lightcoral"), legend=c("Data", "Predicted", "Fitted line"))
qqnorm(resid, col ="cornflowerblue") # Create normal probability plot for residualsqqline(resid, col ="red") # Add a straight line to normal probability plot
💪 Uppgift 4.2
Använd modellen i Uppgift 4.1 för att prediktera koncentration av kväveoxider för observation 10, där employ_dist=10.5857. Beräkna vad residualen blir för denna observation.
Uppgift 4.2 - Svar
Den predikterade koncentrationen av kväveoxider är 0.2662 när employ_dist är 10.5857. Det verkliga värdet är 0.413.
Residualen är 0.413 - 0.2662 = 0.1468
# Write your code herenew_x <-data.frame(employ_dist =c(10.5857))predict(regressionmodel_NOx_employdist, newdata = new_x)
Transformera variablerna i Uppgift 4.1 (avgör själv vilken eller vilka av de två som behöver transformeras). Ett förslag är att använda Tukeys cirkel för att hitta lämpliga transformationer. Anpassa en ny linjär regression med de transformerade variablerna. Utför en modellvalidering (efter transformation) via en residualanalys och kommentera modellens lämplighet jämfört med modellen i Uppgift 4.1.
Uppgift 4.3 - Svar
Vi testar oss fram för att hitta den transformation som ger den rätaste linjen och kommer fram till att vi transformerar både x och y till log(x) och log(y).
Modellvalidering
Jämfört med modellvalideringen i 4.1 ser vi att variansen nu är konstant, vi kan inte se något tydligt mönster i residualplotten. Den är även mer normalfördelad än 4.1 även om residualerna följer samma mönster så är de närmare den röda linjen.
# Write your code hereplot(NOx ~ employ_dist, data = Boston_census_data,col ="cornflowerblue") #otransformerad
# Testar olika transformationerplot(sqrt(NOx) ~sqrt(employ_dist), data = Boston_census_data, col ="cornflowerblue")
plot(sqrt(NOx) ~ employ_dist, data = Boston_census_data, col ="cornflowerblue")
plot(NOx ~sqrt(employ_dist), data = Boston_census_data, col ="cornflowerblue")
plot(log(NOx) ~log(employ_dist), data = Boston_census_data, col ="pink")
plot(NOx ~log(employ_dist), data = Boston_census_data, col ="pink")
plot(log(NOx) ~ employ_dist, data = Boston_census_data, col ="pink")
plot(sqrt(NOx) ~log(employ_dist), data = Boston_census_data, col ="coral")
plot(log(NOx) ~sqrt(employ_dist), data = Boston_census_data, col ="coral")
# Vi väljer att transformera både x och y till log(x) och log(y) för det ger den rätaste linjen.lm_logNOx_vs_logemploy_dist <-lm(log(NOx) ~log(employ_dist), data = Boston_census_data)summary(lm_logNOx_vs_logemploy_dist)
Call:
lm(formula = log(NOx) ~ log(employ_dist), data = Boston_census_data)
Residuals:
Min 1Q Median 3Q Max
-0.19152 -0.07828 -0.01370 0.06108 0.28496
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.218492 0.011524 -18.96 <2e-16 ***
log(employ_dist) -0.327303 0.008832 -37.06 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.103 on 478 degrees of freedom
Multiple R-squared: 0.7418, Adjusted R-squared: 0.7413
F-statistic: 1373 on 1 and 478 DF, p-value: < 2.2e-16
qqnorm(resid, col ="cornflowerblue") # Create normal probability plot for residualsqqline(resid, col ="red") # Add a straight line to normal probability plot
plot(log(NOx) ~log(employ_dist), data = Boston_census_data, col ="cornflowerblue")
log_y_hat <-predict(lm_logNOx_vs_logemploy_dist)orginal_skala_y <-exp(log_y_hat) #transformera tillbaka orginaldata plot(NOx ~ employ_dist, data = Boston_census_data, col ="cornflowerblue")lines(Boston_census_data$employ_dist, orginal_skala_y, type ="p", col ="lightcoral")abline(lm_logNOx_vs_logemploy_dist, col ="lightcoral")legend(x ="topleft", pch =c(0,5, 0,5, NA), lty =c(NA, NA, 0,5), col =c("cornflowerblue", "lightcoral", "lightcoral"), legend=c("Data", "Predicted", "Fitted line"))
💪 Uppgift 4.4
Plotta den anpassade regressionen från 4.3 i icke-transformerad skala tillsammans med observationerna (också i icke-transformerad skala) i en lämplig figur.
Uppgift 4.4 - Svar
Skriv svaret här.
# Write your code herelogy_hat <-predict(lm_logNOx_vs_logemploy_dist) # log scale predictiony_hat <-exp(logy_hat) # original scale predictionhead(logy_hat)
plot(NOx ~ employ_dist, data = Boston_census_data, col ="cornflowerblue") # Data on original scalelines(Boston_census_data$employ_dist, y_hat, type ="p", col ="lightcoral")legend(x ="topleft", pch =c(1, 1), col =c("cornflowerblue", "lightcoral"), legend=c("Data", "Predicted"))
💪 Uppgift 4.5
Använd modellen i Uppgift 4.3 för att prediktera koncentration av kväveoxider i icke-transformerad skala för observation 10, där employ_dist=10.5857. Beräkna vad residualen blir för denna observation. Kommentera resultaten jämfört med Uppgift 4.2.
Tip
Tänk på att ta hänsyn till eventuella transformationer!
Uppgift 4.5 - Svar
Den predikterade koncentrationen av kväveoxider är i icke-transformerad skala ≈ 0.3713.
Residualen är 0.413 - 0.3713 = 0.0417, som kan jämföras med uppgift 4.2 där den var 0.1468. Att residualen är mindre tyder på att modellen passar bättre.
I detta avsnitt ska ni studera multipel linjära regression.
💪 Uppgift 5.1
Anpassa en linjär regression med responsvariabel logaritmerad median_home_value samt förklarande variabler lower_stat_pct och dummy-variabeln borders_charles. Tolka koefficienten för borders_charles.
Uppgift 5.1 - Svar
Från utskriften kan vi utläsa att b_0 (interceptet) är 3.561555, b_1 (lower_stat_pct) är -0.043698 och b_2 (boarders_charles) =0.134560.
Detta betyder att om Borders Charles är lika med 1; alltså att om man gränsar till Charles River ökar det predikterade median_home_value med ≈ 0.135. Om Borders Charles är 0, påverkas inte median_home_value (på grund av 0x0.135 = 0).
# Write your code herelm_median_home_value_vs_lower_stat_pct_borders_charles <-lm(log(median_home_value) ~ lower_stat_pct + borders_charles, data = Boston_census_data)summary(lm_median_home_value_vs_lower_stat_pct_borders_charles)
Call:
lm(formula = log(median_home_value) ~ lower_stat_pct + borders_charles,
data = Boston_census_data)
Residuals:
Min 1Q Median 3Q Max
-0.94794 -0.12591 -0.01925 0.11351 0.89520
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.561555 0.021210 167.916 < 2e-16 ***
lower_stat_pct -0.043698 0.001426 -30.637 < 2e-16 ***
borders_charles 0.134560 0.042267 3.184 0.00155 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2206 on 477 degrees of freedom
Multiple R-squared: 0.6658, Adjusted R-squared: 0.6644
F-statistic: 475.1 on 2 and 477 DF, p-value: < 2.2e-16
💪 Uppgift 5.2
Ni ska nu utforma en modell som predikterar medianhuspriset median_home_value. Ni får endast använda följande förklaringsvariabler:
before_1940_p
crime_rate
radial_access
NOx
dist_fenway_park
Ni får själva välja hur många av variablerna som ska ingå i modellen. Ni får göra vilka transformationer ni vill av variablerna, inklusive responsvariabeln.
Pröva er fram metodiskt när ni väljer vilka variabler ni inkluderar i modellen, och när ni bestämmer vilka eventuella transformationer ni använder.
När ni utvärderar olika modeller kan ni förslagsvis börja med att jämföra adjusted R-squared. När ni med hjälp av adjusted R-squared har identifierat två eller tre modeller som ser lovande ut kan ni utvärdera dessa modeller ytterligare i ett andra steg.
I det andra steget ska ni utvärdera hur väl modellerna predikterar data som inte använts för att anpassa modellen. Ni kan välja en av två alternativa metoder:
Dela in ert dataset i träningsdata och testdata. Anpassa modellen med hjälp av träningsdata, och utvärdera sedan på testdata. Ni kan exempelvis använda de första 350 observationerna som träningsdata och de sista 130 observationerna som testdata.
Använd korsvalidering. Det är en något mer krävande metod, men också något bättre. Ni kan exempelvis göra korsvalidering med 4 folds (4-fold cross validation). Dela då upp ert dataset i fyra delar (del 1: observationer 1-120, del 2: observationer 121-240, del 3: observationer 241-360, del 4: observationer 361-480).
Sortera inte observationerna i Boston_census_data slumpmässigt. Ordningen är redan slumpmässig.
Tip
Tänk på att ta hänsyn till eventuell transformation av responsvariabeln. Om ni exempelvis har valt transformationen \(\log(y)\) är modellens prediktion av responsvariabeln \(\widehat{\log(y)}\). Ni måste då transformera den till \(\hat y\) i responsvariabelns originalskala med formeln \(\hat{y}=\exp\left(\widehat{\log(y)}\right)\). Sedan kan ni räkna ut residualen \(y - \hat y\).
Uppgift 5.2 - Svar
Vi har tränat och testat både model 1 och model 2 med första metod att dela 350 som träningdata och resten 130 är testdata, model 1 har högre RMSE än model 2 därför är model 2 bättre.
# Write your code heredatatrain <- Boston_census_data %>%slice(1:350) # De första 350 observationernadatatest <- Boston_census_data %>%slice(351:480) # De sista 130 observationerna#Använd en linjär regression som modellmodel1 <-lm(log(median_home_value) ~ before_1940_p, data=datatrain)model2 <-lm(log(median_home_value) ~ before_1940_p + NOx, data=datatrain) #Träna modellen 2#Utvärdera model 1y_hatt_test_1 <-predict(model1, newdata=datatest)y_hatt <-exp(y_hatt_test_1)y_test <- datatest$median_home_valueSSE <-sum((y_test - y_hatt)^2)n_test <-130MSE <- (SSE / n_test)RMSE <-sqrt(MSE)RMSE
plot(datatest$median_home_value, resid_model_2, xlab="median_home_value", ylab='Residuals', col ="cornflowerblue")
qqnorm(resid_model_2, col ="cornflowerblue") # Create normal probability plot for residualsqqline(resid_model_2, col ="red") # Add a straight line to normal plot
💪 Uppgift 5.4
Använd modellen i Uppgift 5.2 för att prediktera medianhuspriset för observationerna i datasetet Boston_districts_to_predict (ladda ner).
Uppgift 5.4 - Ladda hem data för prediktion
Ladda in dataseten Boston_districts_to_predict med följande kod.
# Write your code hereload(file =url("https://github.com/StatisticsSU/SDA1/blob/main/assignments/assignment1/Boston_districts_to_predict.RData?raw=true"))
Det här datasetet har endast de förklarande variablerna, dvs inte responsvariabeln. När vi rättar era inlämningsuppgifter kommer vi att jämföra era prediktioner med de faktiska medianpriserna (som vi har tillgång till).
Skriv ut dina prediktioner så att vi enkelt kan se dem när vi rättar.
Tip
Tänk på att ta hänsyn till eventuella transformationer av förklaringsvariablerna!
Totalundersökningen trunkerade medianhusvärdet till 50K för de censusdistrikten som låg över. Vi har tagit bort dessa censusdistrikt. Vi har också tagit bort variabler som är irrelevanta.↩︎
Harrison Jr, D., & Rubinfeld, D. L. (1978). Hedonic housing prices and the demand for clean air. Journal of Environmental Economics and Management, 5(1), 81-102.↩︎
Pace, R. K., & Gilley, O. W. (1997). Using the spatial configuration of the data to improve estimation. The Journal of Real Estate Finance and Economics, 14(3), 333-340.↩︎